﻿using System;
using System.Collections.Generic;
using BarStructures.Deformed;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using BarStructures.Utils;

namespace BarStructures.Analysers
{
    public class FrameAnalyser_XY_DM: Analyser
    {
        private Matrix eigenVectors = null;

        private Vector eigenValues = null;

        /// <summary> Массив смешаных индексов. </summary>
        int[] mergedDegrees;

        public FrameAnalyser_XY_DM(Model model, string name)
            : base(model, name)
        {

        }

        public override void Analyse()
        {
            if (!UseExistingIndexation)
                model.SetAutoSystemIndexes();

            int modelDimension = model.GetNodesCount() * 2;

            //подготовка матриц
            Matrix stiffness = new Matrix(modelDimension, modelDimension);
            Matrix mass = new Matrix(modelDimension, modelDimension);

            #region Построение матриц жесткости и масс
            Matrix elMx, save;
            int barsCount, ind1, ind2;
            barsCount = model.GetBarsCount();
            for (int i = 0; i < barsCount; i++)
            {
                Bar bar = model.GetBar(i);
                ind1 = bar.NodeOne.SystemIndex * 2;
                ind2 = bar.NodeTwo.SystemIndex * 2;

                elMx = ElementHelper.GetGlobalStiffness6DMatrix(bar);
                save = stiffness.GetMatrix(ind1, ind1 + 1, ind1, ind1 + 1);
                save.Add(elMx.GetMatrix(0, 1, 0, 1));
                stiffness.SetMatrix(ind1, ind1 + 1, ind1, ind1 + 1, save);

                save = stiffness.GetMatrix(ind1, ind1 + 1, ind2, ind2 + 1);
                save.Add(elMx.GetMatrix(0, 1, 3, 4));
                stiffness.SetMatrix(ind1, ind1 + 1, ind2, ind2 + 1, save);

                save = stiffness.GetMatrix(ind2, ind2 + 1, ind1, ind1 + 1);
                save.Add(elMx.GetMatrix(3, 4, 0, 1));
                stiffness.SetMatrix(ind2, ind2 + 1, ind1, ind1 + 1, save);

                save = stiffness.GetMatrix(ind2, ind2 + 1, ind2, ind2 + 1);
                save.Add(elMx.GetMatrix(3, 4, 3, 4));
                stiffness.SetMatrix(ind2, ind2 + 1, ind2, ind2 + 1, save);


                elMx = ElementHelper.GetGlobalDescreteMass6DMatrix(bar);
                save = mass.GetMatrix(ind1, ind1 + 1, ind1, ind1 + 1);
                save.Add(elMx.GetMatrix(0, 1, 0, 1));
                mass.SetMatrix(ind1, ind1 + 1, ind1, ind1 + 1, save);

                save = mass.GetMatrix(ind1, ind1 + 1, ind2, ind2 + 1);
                save.Add(elMx.GetMatrix(0, 1, 3, 4));
                mass.SetMatrix(ind1, ind1 + 1, ind2, ind2 + 1, save);

                save = mass.GetMatrix(ind2, ind2 + 1, ind1, ind1 + 1);
                save.Add(elMx.GetMatrix(3, 4, 0, 1));
                mass.SetMatrix(ind2, ind2 + 1, ind1, ind1 + 1, save);

                save = mass.GetMatrix(ind2, ind2 + 1, ind2, ind2 + 1);
                save.Add(elMx.GetMatrix(3, 4, 3, 4));
                mass.SetMatrix(ind2, ind2 + 1, ind2, ind2 + 1, save);
            }

            //DbgPrint.Matrix(mass, "mass_0.txt");
            //DbgPrint.Matrix(stiffness, "stiffness_0.txt"); 
            #endregion

            #region Индексы основных и второстепенных степеней свободы и смешанные индексы

            List<int> primaryDegreesList = new List<int>();
            List<int> secondaryDegreesList = new List<int>();
            for (int i = 0; i < model.GetNodesCount(); i++)
            {
                Node node = model.GetNodeBySystemIndex(i);
                if (node.DegreeSelection.X)
                    primaryDegreesList.Add(i * 2);
                else
                    if (!node.Support.Sx)
                        secondaryDegreesList.Add(i * 2);
                if (node.DegreeSelection.Y)
                    primaryDegreesList.Add(i * 2 + 1);
                else
                    if (!node.Support.Sy)
                        secondaryDegreesList.Add(i * 2 + 1);
            }
            int[] primaryDegrees = primaryDegreesList.ToArray();
            int[] secondaryDegrees = secondaryDegreesList.ToArray();
            mergedDegrees = new int[primaryDegrees.Length + secondaryDegrees.Length];
            Array.Copy(primaryDegrees, mergedDegrees, primaryDegrees.Length);
            Array.Copy(secondaryDegrees, 0, mergedDegrees, primaryDegrees.Length, secondaryDegrees.Length);

            int primaryLength = primaryDegrees.Length;
            int secondaryLength = secondaryDegrees.Length;
            int mergedLength = mergedDegrees.Length;

            #endregion

            //DbgPrint.Matrix(mass.GetMatrix(mergedDegrees, mergedDegrees), "mass_0.txt");
            //DbgPrint.Matrix(stiffness.GetMatrix(mergedDegrees, mergedDegrees), "stiffness_0.txt"); 

            //Рабочие матрицы, без связанных степеней свободы
            Matrix K = stiffness.GetMatrix(mergedDegrees, mergedDegrees);
            Matrix M = mass.GetMatrix(mergedDegrees, mergedDegrees);

            //Вычисление собственных значений и собственных векторов
            EigenvalueDecomposition evd = new EigenvalueDecomposition(M.Inverse().Multiply(K));
            eigenValues = evd.RealEigenvalues;
            eigenVectors = evd.EigenVectors;
        }

        public override DeformedState GetDeformedStateForEigenValue(double eigenValue)
        {
            if (eigenValues == null)
                Analyse();

            List<double> eigenVector = GetFullEigenVecsListForEigenValue(eigenValue);

            if (eigenVector.Count != 0)
            {
                DeformedState deformedState = new FrameDeformedState();
                deformedState.AssignToModel(model);
                int counter = 0;
                double dx, dy;
                for (int j = 0; j < model.GetNodesCount(); j++)
                {
                    Node node = model.GetNodeBySystemIndex(j);
                    dx = eigenVector[counter++];
                    dy = eigenVector[counter++];

                    NodeDisplacement nodeDisplacement = new NodeDisplacement(node, dx, dy, 0);
                    deformedState.Add(j, nodeDisplacement);
                }

                return deformedState;
            }
            else
                return null;
        }

        public override List<double> GetEigenValuesList()
        {
            if (eigenValues == null)
                Analyse();
            List<double> result = new List<double>();
            for (int i = 0; i < eigenValues.Length; i++)
                result.Add(eigenValues[i]);
            return result;
        }

        public override List<double> GetEigenVecsListForEigenValue(double eigenValue)
        {
            if (eigenValues == null)
                Analyse();
            eigenValue = 1;
            int i = 0;
            while (i < eigenValues.Length && !Number.AlmostEqual(eigenValue, eigenValues[i]))
                i++;
            List<double> result = new List<double>();
            if (i < eigenValues.Length)
            {
                for (int j = 0; j < eigenVectors.RowCount; j++)
                    result.Add(eigenVectors[j, i]);
            }
            return result;
        }

        public override List<double> GetFullEigenVecsListForEigenValue(double eigenValue)
        {
            if (eigenValues == null)
                Analyse();

            int i = 0;
            while (i < eigenValues.Length && !Number.AlmostEqual(eigenValue, eigenValues[i]))
                i++;

            double[] evec = new double[model.GetNodesCount() * 2];
            if (i < eigenValues.Length)
            {
                for (int j = 0; j < eigenVectors.RowCount; j++)
                    evec[mergedDegrees[j]] = eigenVectors[j, i];
            }

            return new List<double>(evec);
            
        }
    }
}
